arXiv: 1502.01191v3 [math.DS] 19 May 2015 


Pseudo generators for under-resolved 
molecular dynamics 

What does the Smoluchowski equation tell us about the spatial dynamics of 

molecular systems? 

Andreas Bittracher* 1 , Carsten Hartmann' 2 , Oliver Junge* 1 , and Peter Koltai^ 2 

1 Fakultat fur Mathematik, Technische Universitat Munchen, D-87548 Garching 
2 Institut fur Mathematik, Freie Universitat Berlin, D-14195 Berlin 

May 20, 2015 


Abstract 

Many features of a molecule which are of physical interest (e.g. molecular conforma¬ 
tions, reaction rates) are described in terms of its dynamics in configuration space. This 
article deals with the projection of molecular dynamics in phase space onto configuration 
space. Specifically, we study the situation that the phase space dynamics is governed 
by a stochastic Langevin equation and study its relation with the configurational Srnolu- 
chowski equation in the three different scaling regimes: Firstly, the Smoluchowski equa¬ 
tions in non-Cartesian geometries are derived from the overdamped limit of the Langevin 
equation. Secondly, transfer operator methods are used to describe the metastable be¬ 
haviour of the system at hand, and an explicit small-time asymptotics is derived on which 
the Smoluchowski equation turns out to govern the dynamics of the position coordinate 
(without any assumptions on the damping). By using an adequate reduction technique, 
these considerations are then extended to one-dimensional reaction coordinates. Thirdly, 
we sketch three different approaches to approximate the metastable dynamics based on 
time-local information only. 


1. Introduction 

Langevin and Hamiltonian dynamics constitute established models for the analysis of biomolec- 
ular processes by classical molecular dynamics. While they describe the system at hand 
through the evolution of configuration and momentum coordinates, many properties of inter¬ 
est, such as metastable conformations, conformational transition rates or folding pathways, 
are merely characterized by the configurational dynamics or the dynamics of few collective 
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variables, called reaction coordinates , that span a low-dimensional submanifold of the config¬ 
uration space (see, e.g., mmm)- 

Both from a computational and modeling point of view it is very appealing to describe 
a molecular system just by its position (or reaction) coordinates, since this drastically re¬ 
duces the dimensionality of the problem. Over decades, it has been of major interest to 
derive equations which govern the evolution of these coordinates either exactly [521 [31], or 
approximately with the smallest possible error [II [29]. One popular model for molecular 
dynamics in position space that comes under various names like overdamped Langevin dy¬ 
namics , Brownian dynamics , Kramers equation or Smoluchwski equation is obtained by the 
so-called Smoluchowski-Kramers approximation of the Langevin equation [Ml EH; see also 
emu and the references therein. Yet it is unclear whether there are conditions beyond the 
asymptotic regime of the Kramers-Smoluchowski approximation (i.e. the high-friction limit), 
under which the Smoluchowski equation accurately captures e.g., the folding dynamics of a 
protein in terms of a one-dimensional reaction coordinate. For a general account of this topic 
we refer to [19 ]. 

Aims and scope of this article. In this article, we discuss the accuracy of the Smoluchowski 
equation for the spatial dynamics of a molecular system under various parameter regimes 
where, in each case, our analysis departs from the Langevin equation in phase spaced Our 
presentation of the topic is not claimed to be exhaustive; it rather reflects the authors’ inter¬ 
ests, and their wish to understand how the hierarchies of models used in molecular dynamics 
relate to each other. Parts of this article are based on the recent work [2] by some of the 
authors, however, their analysis in the context of long time scales and in non-Cartesian ge¬ 
ometries is new. The main contribution of this article is that it sketches solutions to answer 
the following questions: 

(a) What is the appropriate generalization of the Smoluchowski equation in generalized (non- 
Cartesian) coordinates, to be used, for example, in reduced-order models of protein folding 
or polymers, and how is it related to a phase space description of a molecular system? 

(b) Is there a closed equation for the spatial dynamics on small time intervals if the underlying 
phase space dynamics is governed by a Langevin equation? 

(c) How well (and in which sense) is a system’s metastable behaviour approximated by the 
Smoluchowski dynamics when the phase space dynamics is generated by Langevin dy¬ 
namics? What are the time scale regimes on which the approximation of the metastable 
dynamics by the Smoluchowski equation may be used? 

The manuscript is organized as follows: Section [2] introduces the basic model of molec¬ 
ular dynamics in terms of deterministic and stochastic differential equations and describes 
an operator-based framework for the evolution of probability densities under these dynam¬ 
ics. This section also introduces the formulations of the stochastic equations in generalized 
coordinates in a non-Euclidean space. Section [3] reviews the concept of metastability based 
on density fluctuations in position space and establishes a connection between Langevin and 
Smoluchowski dynamics on short time scales. A numerically exploitable scheme which replaces 


1 We use the terms spatial, position(al), and configuration (al) interchangeably, when referring to coordinates 
or dynamics. 
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the complicated position space density transport by a rescaled Smoluchowski transport is de¬ 
scribed, along with asymptotic error estimates. Section [4] reviews the approximation quality 
of these methods, gives improved error estimates and discusses the extension to longer times 
scales. A summary and possible future directions are given in Section [5} 


2. Trajectory- and ensemble-based views 

We consider a dynamical system described by d = 3n positional degrees of freedom that 
represent a system of n particles. Let Q cz denote the corresponding configuration space 
and V ( q ) the potential energy of a given particle configuration q e Q where we assume that 
the function V: Q —> M is at least twice continuously differentiable, polynomially growing at 
infinity and bounded from below. 


2.1. Models for molecular dynamics 


We introduce three typical models for molecular dynamics. The simplest model to describe 
the motion (qt)t^o of the particles in vacuum, i.e. without external influences like a solvent is 
given by Hamilton’s equations 


ft -~V a H( q ,p), 


( 1 ) 


where p e V = denotes the vector of conjugate particle momenta, and 


H = -p-M ] p + V(q) 


is the Hamiltonian (total energy) of the system, with M = diag(mi,..., md) denoting the 
mass matrix. Depending on the type of system or when transformed to generalized coordi¬ 
nates, the mass matrix M can be a general symmetric positive-definite, possibly position- 
dependent matrix. 

In the presence of a heat bath or solvent, one typically adds a drift-diffusion term, arriving 
at the Langevin equation , 

j- - V p H(q,p) 

1 (2) 

— = - X7 q H(q,p ) - 7 X7 p H(q,p) + a( t . 

The term — 'y\7 p H = —'yM~ 1 p, with 7 6 M. dxd being symmetric positive definite, models 
the drag through the solvent, the term aQ accounts for random collisions with the solvent 
particles [35], Chap. 4], Here, 0 is an uncorrelated, zero-mean white noise process that can 
be formally interpreted as the (generalized) derivative of a standard d-dimensional Brownian 
motion, and aa T e M dxd is the noise covariance matrix. In order to keep the system at a 
constant average kinetic energy, damping and excitation have to be balanced, which is ensured 
by assuming that noise and drag coefficients satisfy the fluctuation-dissipation relation 


27 = /3<Tcr r , 
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where f3 > 0 is the inverse temperature in the system. Choosing 7 or a is a modelling issue 
and thus depends on the particular problem at hand. As we will see later on, both 7 and a 
may even be position dependent. 

If the friction in the system is uniformly large, i.e. v ■ r yv » v ■ Mv for all v e M rf , the 
Langevin equation can be replaced by the Smoluchowski equation 

f§ = -W ( <?) + ^, (3) 

or, using the common notation for Ito stochastic differential equations (see m ), 

7 dqt = -VI f(q t )dt + odw t , 

where wt is a standard Brownian motion in The Smoluchowski equation is also termed 

overdamped Langevin equation and can be derived from ([2]) by letting v ■ 'yv —> 00 under a 
suitable rescaling of time; for a precise statement, see 135, Theorem 10.1] or the derivation 
given below on pp. [5]-[6j 


In some cases a description of the stochastic dynamics in a different coordinate system 
is needed. The aim of this subsection therefore is to derive the Smoluchowski equation in 
generalized coordinates. This is most conveniently done by resorting to the canonical form of 
the Langevin equation ([2]) that we will state first. 


Langevin equation in generalized coordinates. To state the Langevin equation ([2]) in canon¬ 
ical form, we consider a diffeomorphism d>: Q' —> Q between configuration spaces that has a 
cotangent lift T*<b: X —* X’ given by 

(?,p) *-*• ($ _1 (?), ((Vd>o $" 1 )( g )) T p) . 

Using Ito’s formula [ 361 Theorem 4.2.1], the Langevin equation ([2]) can be written in the new 
configuration variables u = < k~ 1 (g) and their conjugate momenta v = ((V<£ o 4>” 1 * * )(g)) T p: 
introducing the new (possibly position-dependent) drag and noise coefficients by 

7 = V<f> T 7V<I>, (T = V4> T cr, (4) 


the Langevin equation can be recast as mm 

du „ . 

* = (<■>") 


dv 

dt 


(5) 


= -X7 U H(u, v) - 7 (u)V v H(u, v) + a(u)Ct ■ 


Here H denotes the push-forward of the Hamiltonian H to the new coordinate system, 

H(u,v) = X -v ■ (G(u)) _1 u + V{u ), 


with V = V o (E^ 1 and G = VT 7 AfV4> being the mass metric tensor induced by the trans¬ 
formation <h. 

It can be readily seen that, when the original drag and noise coefficients satisfy the 
fluctuation-dissipation relation, then so do the transformed coefficients: 


27 = (3dd T . (6) 

2 In the following, we will interpret stochastic differential equations such as ([2]) or ^ in the sense of Ito, 

which implies that, for stochastic differentials such as dqt, a generalized chain rule known as ltd’s formula 

or Ito’s Lemma applies HD Theorem 4.2.1], 
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Derivation of the Smoluchowski equation in generalized coordinates. It is now possible 
to derive the Smoluchowski equation in generalized coordinates from the canonical Langevin 
dynamics ([ 5 j using formal asymptotics. To this end, let us scale the original drag and noise 
coefficients according to 7 >—> 7/e and <r >—> cr/yfe where £ > 0 is a small parameter. Clearly, 
the scaling preserves the fluctuation-dissipation relation <©> and it leads to the Langevin 
equation 


— = -V u H(u, v) - 7 (u)V v H(u, v) + —a(u)Ct ■ 

at e ye 


( 7 ) 


For notational convenience, we have dropped the twiddle signs on the transformed Hamilto¬ 
nian H and the coefficients 7 and a. To study the e —* 0 limit of ([T]) we seek a perturbative 
expansion of the associated backward Kolmogorov equation^ 


d t (/) e (u, v, t ) = A Lan (j) £ (u , v, t ), <fi £ (u, v, 0 ) = v ) 


( 8 ) 


following the approach described in [371 [ 39 ]. To begin with, we notice that the backward 
operator A^ an in ([8]) admits the decomposition (see also p. [ 9 ] below) 

^Lan = ^Ham + “^OU , 

£ 


with 


^Ham = y v H ■ V u - V U H • V„ 


and 

-4ou = \ a(jT '- - in G ~ lv ) ■ 

We consider a perturbative solution of ([8]) that is of the form 

<f) £ = (f)Q £(f> 1 + £ 2 </>2 + . . . 


with cf>i = 4>i(u, v , t). Inserting the ansatz into the backward equation and equating powers of 
e we obtain a hierarchy of equations, the first three of which read 


^OU^O 

= 0 

( 9 ) 

^OU^l 

= dt 4 > 0 — Hnam^O 

(10) 

^OU<^2 

= dt 4 > 1 — • 

(11) 


Note that Hou is a second-order differential operator in v with u appearing only as a pa¬ 
rameter. By the assumption that 7(-) is symmetric positive definite with uniformly bounded 


3 The Kolmogorov backward equation is a partial differential equation governing the evolution of observables. 
Specifically, if (Xt)t.^o is the solution of an Ito stochastic differential equation, such as 0- then, for any 
integrable function (f>o: X —* R, 

=E x [MXt)] 


satisfies = A<j> with initial condition 0) = (f>o{x). Here A is the infinitesimal generator associated 
with the process (W)t^o, and E^f ] denotes the conditional expectation over all realizations of the process 
starting at Xo — x. We introduce the dual concept of forward equation and the corresponding generator 
cf. also Section 7.3]. 


in Section 


2.2 
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inverse, implies that 0o does not depend on u. By a standard closure argument (a.k.a. cen¬ 
tering condition ), it thus follows that dt<j >o = 0. 

Closely inspecting the resulting equations the next nontrivial term, cj) i, is found 

to be the solution of the backward equation 

dt'l/j = I (^dHam 7 4oU^HamV ; ) Qu{v)dv , (12) 

JR d 


where tp = 4>(u,t) is independent of v, and g u is the solution to Aq^q u = 0, with Aq^j 


being the formal L'? adjoint of Aou- Equation (12) must be equipped with a suitable initial 
condition 0) = tpo(u). 


Before we evaluate the right hand side of (12), a few remarks are in order: 


1. The function g u in (12) is the unique invariant probability density with respect to the 


Lebesgue measure on the momentum space V = T* Q' (that we can identify with 
of the Ornstein-Uhlenbeck process generated by Aqu ■ It is given by 


8u(v) 


7 ^) (det(G(u))) 1/2 exp ^-^v(G(u)) T v) 


(13) 


and satisfies AQ V g u = 0. 


2. The inverse of the operator Aqu in (12) is only unambiguously defined when it is acting 
on functions that are in the range of Aou- By the Fredholm alternative (see m 
Sections 3.12 &; 4.3]), the range of a linear continuous operator on some Banach space 
is the orthogonal complement of the kernel of its adjoint: 


ran(A 0 u) = (ker(Aou)) = V o , 


where 


V 0 :=\feC 


(V): j f(v,u)eu(v)dv = o| c L\{X') 


That is to say that the linear equation Aoud* = c has a solution, if and only if c averages 
to zero under g u . Note that c = ^HamV’ is linear in the momenta v , hence A HamV 7 6 E). 


As a consequence, A q Jj in (12) is well defined. 


The formal expansion ([9])-( 12) now suggests that the solution of the Langevin-based backward 
equation (l8|) and the solution to the limiting system (12) satisfy 


'(V) -ip(-,t/e)\\ v —* 0 , 


0. 


(14) 


for some suitable norm on V cz Indeed, standard results from homogenization theory 
for parabolic partial differential equations (e.g. |37 ( 138 ] ) imply that, under certain regularity 
assumptions on the coefficients, the convergence is uniform in X 1 x (0, T ) for any finite T 
(cf. also |35j Theorem 10.1]). 


As we show in the appendix, the operator on the right hand side of (12) reads 

A = /8 _1 A -W-V, 


(15) 
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where 


1 


V = 7 1 V and A = 


V det 7 


V 


(Vdet77 , 


denote gradient and Laplace-Beltrami operator with respect to 7 . The differential operator 
A has a straighforward interpretation as the infinitesimal (backward) generator of the Srnolu- 
chowski dynamics on the configuration space Q! c X', with the position dependent drag 
matrix acting as metric tensor on Q!. Alternatively, one may regard A as the generator of 
the Smoluchowski dynamics on a Riemannian manifold Q] endowed with the metric tensor 
h = and a position dependent friction matrix 7 . Our findings are summarized in the 

next lemma. 


Lemma 1. The Smoluchowski equation in generalized coordinates reads 


du 

l(u)—=-VV(u) + g(u) + a(u)( t , u 0 = u , 


where g = (< 71 ,... ,g n ) has the entries 

1 1 d 

~R Z J . /Uh 


9i = 


P ji ' Vdet 7 du k 


(y/defy 7 ^ . 


(16) 


(17) 


Remark 1. The additional drift term g in the Smoluchowski equation is due to the geometry 
of the configuration manifold Q! and the interpretation of the Smoluchowski equation in the 
sense of ltd. Formally, it can be seen to be related with the first order derivative in the 


expression for the Laplace Beltrami operator in (15). 


Remark 2. The Smoluchowski equation (16) in generalized coordinates likewise follows by 
transforming the original Smoluchowski equation & in Cartesian coordinates to the new 
coordinate system using ltd’s Lemma (36. Theorem 4-2.1]. As a consequence, the stochastic 
convergence of the spatial component of the high friction Langevin equation to the solution of 
the (time-rescaled) Smoluchowski equation that has been proved in 135k Theorem 10.1] should 
be inherited by its non-Cartesian counterpart. 


2.2. The transfer operator 

We shall now examine how probability densities evolve under the Langevin or the Smolu¬ 
chowski dynamics. To this end, call x = (q,p) or x = q the state vector, depending on which 
type of dynamics is considered, and X = Q x V or X = Q the phase space or state space, 
respectively. For any given xq e X, we seek the probability density ft of xt for some t > 0 
with respect to the natural (Liouville or Lebesgue) measure on X. Now let B c X be any 
measurable subset of X, and define the stochastic transition functiorj^Jof the dynamics (xt)t^o 
by 

p(t , x, B) = Prob \xt e B \ xq = x\ , (18) 


Further call 


Pn(t, A, B) \= Prob M [x t e B \ x 0 e A] 


(19) 


4 It is common to denote the transition function and transition probabilities by p. We hope that the clash in 
the notation with the conjugate momenta is not going to confuse the reader, since the transition function 
and -probabilities are always going to be functions of three variables. 
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the transition probability between the measurable sets A a X and B cz A, where Prob^ 
indicates that the initial condition is distributed according to a probability measure p. 
For the long term macroscopic behaviour of the system, sets A play an important role 
for which p^(t, A, A) « 1 for some physically relevant measure p and some characteristic 
times t > 0. 

The transition probability p(t,x,B ) can be derived from a stochastic transition kernel or 
transition density T via 

p(t, x, B) = f <l>(t,x,y)dy, (20) 

Jb 


where T is the fundamental solution of the Fokker-Planck equation (26). Existence and 
uniqueness of the transition density T follow relatively under mild conditions (see e.g. |28L 
Chap. 11.7]). 

Now let an initial density /o = dp/dx be given. The density ft describing the distribution 
of the system at time t > 0 is then given by 


ft{y) = f Mx)^(t,x,y)dp(a 

Jx 


( 21 ) 


Equation (21) can be seen as the definition of the transfer operator with lag time t: 

P l f o(®) := ft(x). 


By linearity we can extend the definition of P t from probability densities to arbitrary inte- 
grable functions, and it will be convenient in what follows to consider the transfer operator 
as a family of linear maps P l : ( X ) —> C l ^{X). This family of linear operators has the 

Chapman-Kolmogorov (or semigroup ) property: 

(i) P°f = /, 


(ii) P t+S f = P\P s f) for all s,t > 0. 

We also have non-expansiveness in the induced operator norm, ||.P < || ^ 1, and positivity, 
P f f ^ 0 for / Ss 0. 

The transition probabilities (19) can be conveniently expressed in terms of the transfer 
operator. If we define the scalar product on the space C 2 ^{X) of square integrable functions 
by 

(f,g)u= f f( x )g( x )dp(x), 

Jx 

then, for any measurable set A with p(A) > 0, 


Pn(t, A, B) 


1 

p(A) 


L 


P l XA dp 


1 

p(A) 


L 


P t XAXB dp 


<■ pt XA1 XB)ti 

{XA, XA)n 


with % being the indicator function. 


The forward generator. The semigroup property means that P < is “memoryless”, i.e. that 
0 @ generate a Markov process. Noting that 

P t = f^pt/ny 
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we may conclude all relevant information about the density transport is already contained in 
P T for arbitrarily small r. This is formalized by looking at the operator 

P T f - f 

Lf = lim - — (22) 

t—> o r 

that is defined for all /, for which the limit exists. L is called the forward generator or 
infinitesimal generator of P t . 

From its general form for arbitrary Ito diffusions [26j p. 282], we can derive the generator 
for our dynamics. For the Hamiltonian dynamics ([Tj) and functions / 6 where C 1 is 

equipped with the supremum norm, the operator L is given by 

^Ham = Vqiif ■ V p — V pH ■ V q , (23) 

where the dot denotes the Euclidean inner product, and V 9 , V p are the gradients with respect 
to q or p. In case of Langevin dynamics ([2]) and / e C 2 (X), we have 

^Lan = ^Ham + \ a(jT : V p + 7 ^pH • V p + 7 : , (24) 

where the notation A: B := tr (A T B) denotes the inner product between matrices A,Be 
R dx , and \/ 2 denotes the matrix of second derivatives with respect to p. Finally, the generator 
of the Smoluchowski dynamics © reads 

^Smoi = /TV 1 : V 2 + ( 7 _1 V,K) • V 9 + 7- 1 : V 2 V , (25) 

with S7 2 being the matrix of second derivatives in q. and we have used that 27 = (3acr T . 


Fokker-Planck equations and invariant measures. By definition of the forward generator, 
the evolution of probability densities ft associated with any of the stochastic dynamics © © 
is described by a parabolic transport equation of the form 


St ft = Lf t , f t=0 {x) = g(x) , 


(26) 


that are called Kolmogorov forward equations or Fokker-Plack equations [26|, with L being 
either L^ an or -Lsmol- When 7 = 0, then the Fokker-Planck equation with L^ an turns into the 
Liouville equation that describes the transport of probability densities under the Hamiltonian 
dynamics ([Tj). 

Probability measures that are invariant under the dynamics play a prominent role. The 
corresponding densities are fixed points of P f for any t ^ 0, and equation (22) implies that 
they lie in the kernel of L. For the stochastic processes considered here, the invariant density 
can be shown to be unique (cf. [33J)- For the Langevin dynamics ([2]), the unique invariant 
probability density is the canonical density 


fam(q,p) = -^exp(-/3 H(q,p)) 

= ^ ex P ^~ 2 P ' M ~ lp 


=-fv(p) 


^exp {-pv(q)), 

v j 


(27) 
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with Z = Z V Z Q and 


Z V = exp (^-^P ■ M l p^ dp, Z Q = l exp ( - f3V(q))dq . 

For the Smoluchowski dynamics (§, the unique invariant measure has the density fQ(q), 
which is called the Boltzmann density or Gibbs-Boltzmann density. (We assume throughout 
that exp(— f3V) is integrable). 

Under fairly mild assumptions on the potential V, the invariant densities can be shown to 
be the unique asymptotically stable fixed point of P t , which implies that P t fo converges to 
the stationary distribution for any initial density /o [31] • The Liouville equation associated 
with the Hamiltonian dynamics 0 is known to have infinitely many stationary solutions, 
among which is f can . 


Remark 3. It can be readily seen that the Smoluchowski dynamics (16) in generalized coor¬ 
dinates us Q! has the unique invariant probability measure 


dp(u) = (f Q o $)(u) dE(u), 


(28) 


with 


dTi{u) = \J det h{u) du 

being the Riemannian volume element on Q! where h(u ) = (V < f >:jr V<l ) )(u) is the corresponding 
metric tensor. Note that (28) is simply the pullback of the Boltzmann distribution in Cartesian 
coordinates by the coordinate transformation <f>. As a consequence, dp/du is the Q!-marginal 
of the canonical density f can ■ To see this, replace the metric tensor h on Qf by the generalized 
mass matrix G = V$ T MV$ or the corresponding expression for the friction coefficient 7. 
This does not change the invariant measure as the constant mass or drag matrices cancel out. 


2.3. More on semigroups and their generators 

Before we proceed, let us recall two results relating the transfer operator semigroup and its 
generator. For our purposes the main connection between them is given by the following 

Theorem 1 (Spectral mapping theorem [3D]). Let X be a Banach space, T l : X —> X, t ^ 0, 
a Cq semigroup of bounded linear operators (i.e. T 4 / —* f as t — > 0 for every f e X , and T l 
bounded for every t), and let A be its infinitesimal generator. Then 

e ta.{A) c ^.( 7 *) C e ta ^ A) u {0}, 


with a, denoting the point spectrum. The corresponding eigenvectors are identical. 

Evidently, a function / is an invariant density of P l for all t ^ 0, if and only if Lf = 0. 
Further, since ||P*|| 1, the eigenvalues of L lie in the left complex half-plane. The family 

P f can be approximated by a truncated “Taylor series”: 

Example 1 ([2]). If f is 2 N + 2 times continuously differentiable and L n f, n = 0,1,..., N, 
is square-integrable with respect to p, then 


N 


p ‘f - 2 -pf c2 - °( (iV+1 ) f° r f - 0 . 


71=0 


where C 2 denotes the C?-norm with respect to p. 
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3. Spatial dynamics and metastability 

Consider an infinite number of systems modeled by ([2]) in thermodynamic equilibrium, i.e. 
identically and independently distributed according to / can (this collection of systems is called 
an ensemble in statistical mechanics). We are now interested in the portion of these systems 
which undergo a certain configurational change, i.e. leave a subset A a Q and enter another 
subset B cz Q. For this, we track the evolution of XAfcan, which is given by P t (xAfc an)- This 
will be the starting point of the subsequent analysis. 

The spatial transfer operator. Since we are only interested in the distribution of their 
configurations at time t ^ 0, we compute the marginal with respect to q. The resulting 
spatial transfer operator for some u e £?j a is @2), :.50j 

S tu ( q ) : = ^ p Ln{u(q)fcan(q,p)) dp. (29) 

Metastability on configuration space. Using the scalar product 

(u,v)f Q :=\ u(q)v(q) f Q (q)dq 
JQ 

(which gives rise to the norm || • ), and the “slice” T(A) := {( q,p ) 6 Q \ q e A} in state 

fQ 

space, we define transition probabilities between slices via 

V (t, T(A),T{B)) = XA ’ X ^ >/a ■ (30) 

\XAi XA/f Q 

We call a disjoint union Ai u ... u A n = Q of position space metastable or almost invariant 
if 

p(t,T(A j ),T(A j )) » 1, j = l,...,n 

for the time scales t > 0 of interest. Other, more sophisticated, notions of metastability can 
be found in pf3j| . 

The link between almost invariant/metastable sets and eigenvalues close to one and the 
corresponding eigenvectors of some transfer operator was first established in [8] and used 
for conformation dynamics in [9]. We here cite an extension to a broader class of transfer 
operators from [251 . 

Theorem 2 (Application of [25], Theorem 2). Let a(S t ) c: [a, 1] with a > — 1 and X n ^ ^ A2 < Ai = 1 

be the n largest eigenvalues of S f , with eigenvectors v n ,. .. ,v\. Let {Ai, ..., A n } be a mea¬ 
surable partition of Q and II be the orthogonal projection onto span(yyin ■ • ■ iXA n )- Then 

1 + P 2 X 2 + ■ ■ ■ + p n Xn + c ^ p(t, T(Ai), T(Ai)) + • • • + p(t, r(A n ),r(A n )) ^ 1 + A2 + • • • + A n , 

where pj = ||IIujj| e [0, 1] and c = o( 1 — p 2 + ■ ■ ■ + 1 — p n )• 

Remark 4. We briefly discuss some specific properties of the spatial dynamics that are useful 
to understand the concept of pseudo generators outlined below. 
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The spatial transfer operator S'* from (29) satisfies the assumptions in Theoremsee n 
Appendix B]. Unfortunately, S t lacks the semi-group property and so cannot he the 
solution operator of an autonomous transport equation like the Fokker-Planck equation. 
Equivalently, spatial dynamics is not induced by an ltd diffusion process and thus has 
no infinitesimal generator in the sense of (26). 


2. The closer the eigenvalues A2, ■ ■ ■, A n are to 1, the more metastable can a partition po¬ 
tentially be (upper bound in the theorem). How metastable a given partition is, is con¬ 
trolled by the pi (lower bound), which measure the constancy of the eigenfunctions on 
the partition elements. The better the eigenfunctions can be approximated by piecewise 
constant functions over the partition, the closer the pi are to 1, and the more metasta¬ 
bility is guaranteed by the lower bound. Also, note that since S* is not a semigroup, the 
eigenfunctions Vi depend on t (cf. Figure [7] below). As a consequence, metastability of 
a partition must be understood with respect to the characteristic time scale t > 0; 


3. It is not necessary that the sets Ai form a full partition, i.e. UA = Q- Similar results 
to the above have been obtained for non-complete partitions, where the Ai are considered 
to be cores of the metastable sets; cf. JJdj Section 5]. 


3.1. Pseudo generators 


Even though the spatial dynamics is lacking the semigroup property, it is possible—at least 
formally and in analogy with (22)—to differentiate S t at t = 0. We will see in the following 
that the resulting operators can play the role of the infinitesimal generator in the context of 
metastability analysis. 


Definition 1. Let A! be a Banach space, T* : X —► X, t > 0 be a time-parametrized family 
of bounded linear operators. The operator 


d_ 

dt 


T t f = 


lim 


'j-'t+h j? _ j~*t jt 


H§^ Tt ) th en-th 


For T* = P f , the transfer operator of an Ito process, the pseudo generators are simply 
G n = L n , where L is the infinitesimal (forward) generator. 

The pseudo generators of the spatial transfer operator S* can be expressed by the generator 
Than of the full Langevin transfer operator: 

Lemma 2 (pj). The n-th pseudo generator G n of S t takes the form 

G n u(q ) = ml L ' LMq>f can ( q,p )) dp- 

Explicitly, we have 
(1) G 1 = 0, 


is called the time-derivative of T*. Iteratively, we define by ^T* : = 
time-derivative. The operator 


r ■= — _ r t \ 

n ■ dt nT 


t =0 


is called the n-th pseudo generator o/T*. 
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(2) G*2 = —A — VU ■ V. In particular, G 2 is independent of 7. 

Surprisingly, one has 

Corollary 1 ([2]). The pseudogenerator G 2 (of the spatial transfer operator) is the infinites¬ 
imal generator of the Smoluchowski dynamics: 

G 2 = Gsmol- 


Remark 5. Note that G 2 = Gg mo i has the form of the backward Smoluchowski genera¬ 
tor Ag mo \ (cf. Section^). Still, G 2 is also the forward generator of the Smoluchowski process, 
if distributions are thought of as distributions with respect to the Gibbs-Boltzmann density fg. 
This is in accordance with the definition of the spatial transfer operator (29), which also de¬ 
scribes redistribution of mass with respect to fg. The formal coincidence “Gg mo i = Ag mo i ” is 
not accidentally, but rather it reflects the reversibility of the Smoluchowski process. 


Taylor reconstruction of the spatial transfer operator. It is natural to ask whether there is 
an analogue of Proposition [I] for S' 4 and its pseudo generators. We have the following result: 

Theorem 3 (j2j). If u is sufficiently regular, then 


K fk 


s ‘“ - S 


k =0 


= 0{t K+l ), (t —0). 


£2 

fQ 


Unfortunately, for k > 3, higher derivatives of the potential V appear in the expressions for 
Gk, which are therefore impractical to work with, while the gradient VU is typically available. 
We call 


t 2 n 


R t u := ^ id +—G 2 ^u = u + — A u — Vu ■ VU 

the (2nd order) Taylor approximation of S t such that if u is sufficiently regular, 

ILs^u - i^ttlLa =o(t 3 ), (t —0). 

L fQ 


(31) 


Exponential reconstruction. Unfortunately, unlike S f , R l is not norm-preserving for den¬ 
sities with respect to / can . Therefore, when transporting u, we lose the interpretation of 
(Bfu) / C an as a physical density. Moreover, R l u is not even bounded in t [2j. This quickly 
(i.e. for small t) destroys the interpretation of the eigenvalues of Rf as a measure of metasta¬ 
bility. 

An alternative approximation to S t preserves those properties. The Taylor approxima¬ 
tion (31) already suggests that 5* behaves similarly to a ^t 2 -scaled Smoluchowski dynamics. 
Hence, we define 

E‘f C‘l 2 ,/. 


where Pg mol is the semigroup of transfer operators generated by G2. This operator is integral¬ 
preserving with respect to the weight fg [2], and we get the following analogue for Proposi¬ 
tion 111 
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Lemma 3 (|2, Lemma 4.10]). If u is sufficiently regular, then for t —> 0, 


In particular, 



N 


L 

n= 0 


n! 


u 



0(t 2N+1 ). 


E t u-S t u La 
L fa 


0(t 3 ) (t 


0 ). 


Reconstruction of eigenspaces. The error asymptotics carries over to the spectrum and 
eigenvectors of S t , R l and E l in the following way: 

Corollary 2 ([2]). Let u be a sufficiently regular eigenvector of R l or of E l to eigenvalue A. 
Then 

115*1* - AuL.2 = 0{t 3 ). 
fa 

Thus, for small t we may interpret dominant eigenpairs (u, A) of E t and R t as good ap¬ 
proximations to dominant eigenpairs of 5*. Hence, they can be used to define metastable sets 
following the spatial decomposition approach in m- The eigenfunctions of interest, those 
of S t , E t , and i? 4 , can be shown to be sufficiently regular under fairly general conditions, 
cf. pi Appendix C]. 

Remark 6. Corollary is in accordance with functional analytical results by Nier and co¬ 
workers (e.g. \2ff ) that show that the dominant spectrum of the non-reversible Langevin dy¬ 
namics is real-valued and close to the spectrum of the reversible Smoluchowski dynamics, 
even for moderate values of the friction coefficient. In \24f , this somewhat surprising result is 
obtained by large deviations arguments for the small-noise limit using the Witten Laplacian 
representation of the (hypoelliptic) Langevin generator, whereas the considerations here and 
in are based on small-time asymptotics of the spatial Langevin dynamics. We believe that 
a connection between these results is that the small-noise limit can be understood as an expo¬ 
nential rescaling of time as is suggested by large deviations theory; cf. m We refrain from 
going into details here and leave the analysis of this interesting connection to future work. 


3.2. Towards spatial generators in essential coordinates 

We have discussed the concept of the spatial transfer operator that is obtained by projecting 
the phase space dynamics onto the spatial components. We shall now consider the restriction 
of the dynamics to a given collective variable, also termed essential coordinate. To this end, 
let £: Q —> Z c M be a smooth map with the property that, for every regular value z e Z of 
£, the level sets 

M z = {q 6 Q: £(q) = z} cz Q 

are smooth submanifolds of Q with codimension 1 (i.e. hypersurfaces). We suppose that £ is a 
physically relevant observable of the dynamics, such as a reaction coordinate or some collective 
variable that monitors a conformational transition, and call £ the essential coordinate ; the 
unessential coordinates are then implicitly defined through the foliation of Q by the map £, 
in other words: the unessential coordinates parameterize the leaves A4 Z of the foliation for 
every (regular) value z of £. 
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To define the analogue of the spatial transfer operator (29) for the essential coordinate, 
firstly note that HQ Section 3.2] 


f g(q)dq= f (\ 5 |V^| 1 da Adz 

JQ Jz \Jm z / 


(32) 


for any integrable function g: Q 


where da, denotes the Riemannian volume element on 


A4 Z . Equation (32) is called the coarea formula and can be considered a nonlinear variant of 
Fubini’s theorem. 

Together with the law of total expectation, the coarea formula thus entails that the canon¬ 
ical probability measure p conditional on £(q) = z has the form 


dp z = T7rr/can|V£| 1 da z dp , 

N(z) 


with the normalization constant 


N(z) = f /can|V^| 1 da z dp. 
JM,xV 


The spatial transfer operator S 4 for essential coordinates can now be defined as 

SlsM Z ) : = f -PLan( u; (?(9))/can(g,P))|V^(g)|' 1 (io- x: dp, 

iv \ z ) JMzxV 


(33) 


(34) 


(35) 


Projected pseudo-generators. To compute the corresponding pseudo-generators, let p be 
the configurational marginal probability measure that is obtained by projecting p onto the 
configurations by integrating out the momenta, i.e., dp(q) = fQ{q)dq. Let us further introduce 
a projection operator U z : C 2 p (Q) —> £ 2 (Q) by 

( n Z u)(z) = 1 f u(q) faiq^Ciq^da^q) (36) 

^Q\ z ) JQ 

where /q is the g-nrarginal of / can and Nq is the corresponding normalization constant for the 
conditional density. It can be readily seen that, II- is an orthogonal projection with respect 
to the natural (weighted) scalar product in the space C 2 p {Q) and amounts to the expectation 
of functions with respect to p conditional on £(q) = z. 

Thus, for functions u(q) = w(£(q)), the reduced spatial transfer operator 5 ass and the 
spatial transfer operator p9] ) are related by (cf. ]42j) 

S lss w ( z ) = (^ Z s t u)(z) . (37) 


The last identity is helpful in computing the corresponding pseudo generators G® ss . Here 
we are interested only in the second pseudo generator G| ss , for which we have the following 
analogue of Lemma [2j 


Lemma 4. For sufficiently smooth functions u(q) = w(£(q)), the n-th pseudo generator of 
5g SS reads 

G e f s w(z) = ( U z G n u)(z ) 


Specifically, we have 


ar-r'nz^+Hz)^. 
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with the noise and drift coefficients 

a{z) = (IE|V£| 2 )(z), b(z) = (iM/r'A? - v£ • W)) (z ). 

Proof. The first part of the assertion is a straight consequence of Lemma [2] and the coarea 
formula. As for the second part, observe that the second pseudo generator is given by G 2 = 
/3 _1 A — W ■ V which by chain rule implies: 

G 2 io(£(g)) = r'lVei 2 w"{z)\ z= ^ q) + tf-'AZ - V£ • VV)w\z )\ z=aq) . 

Letting the projection 11^ act from the left using that U. z w'(z)\ z= ^ = w'(z ) and likewise 
n * w”(z) |. = ^( 9 ) = w"(z ) gives the desired result. □ 

A few remarks are in order: 


1. In accordance with Corollary [lj the second projected pseudo generator G| ss is the in¬ 
finitesimal generator of the diffusion 


dz 

dt 


= h(z) + V 2 P t , 


(38) 


with a(z) = yafz) and Q being a one-dimensional uncorrelated Gaussian white noise 
process. Equation (38) has been derived by Legoll and Lelievre [29] using first-order 


(Markovian) optimal prediction. 


2. In 


the authors prove an error bound for the projected dynamics (38) under the 


assumption that the conditional probability fj, z satisfies a logarithmic Sobolev inequal¬ 
ity. We refrain from transferring the analysis to our situation as logarithmic Sobolev 
constants are difficult to estimate (beyond the case of strictly convex potentials or in 
the zero-noise limit), hence the approach is of limited practical use. Nonetheless, we 
believe that the projected pseudo generator G| ss will provide a good approximation of 
the dominant spectrum of LLan whenever £ is a slow coordinate relative to the unessen¬ 
tial configuration variables and the momenta. 


3. If |V£| is bounded above and away from zero, it can be shown (see [29]) that the process 
{ z t)t:no generated by L ess = G| ss has the unique invariant measure 

dv(z) = exp (—(3F(z))dz 

with 

F(z) = -/3- 1 logjf Q \VZ\- 1 da z 

being the thermodynamic free energy in the essential coordinate. Note that u = /ro£ _1 
is the push-forward of the canonical distribution by £ (i.e. the £-marginal). Naively, one 
might expect the projected Smoluchwski equation to be of the form 

^ = -F\y) + y/2 , (39) 
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and it can be shown that (38) can be transformed into (39) according to y = cp(z) using 
Ito’s Lemma with cp being the volatility transform 


<p(z) = f (er(s)) 1 ds 

Jo 

that leads to a Smoluchowski equation with unit noise coefficient 
quence, (38) can be equivalently expressed as 

^ = -a(z)F'(z) + ^~ 1 a , (z) + V 2/3~ 1 cr(z)(t , 

which is exactly the one-dimensional analogue of (|l6[) with 7 = a 1 . 


As a conse- 


(40) 


4. In order to use G| ss in metastability analysis (analogous to G 2 in section 3.1), it has 
to be discretized. The method of choice is spectral collocation due to the regularity of 
the objects of interest [2] (i.e. eigenfunctions of Sg ss ). Here, collocation requires the 
evaluation of G^ s (pi{zj) for ansatz functions (pi at collocations points Zj (see Section 
[5] for details). This in turn requires the evaluation of the noise and drift-coefficients 
a(zj), b(zj) in Lemma [4j which involve (potentially high-dimensional) integrals that 
represent averages over the non-essential degrees of freedom; see, e.g., [HI \22l [HO] for 
Monte-Carlo methods to efficiently compute these high-dimensional integrals. 


5. We should stress that Lemma [4] can be readily generalized to multidimensional reaction 
coordinates, however, in general (expect for the case of pairwise orthogonal reaction 
coordinates) it is unclear whether the physical interpretation of the projected equation 
as a reversible diffusion in the free energy landscape is retained. 


4. Approximation quality for larger time scales 


^ 2 12 , 

We have seen in Section 3.1 that E l = P<A o] approximates S* well (pointwise) for small 
times t. However, for metastability analysis, spectral properties of the spatial operator for 
larger time scales are of interest. In this section we make use of two well-known techniques— 
perturbation expansion, already seen in Section 2.1, and the Mori-Zwanzig formalism —with 
the aim of explaining the approximation quality of pseudo generator reconstructions of S l , 
and extending them to larger time scales. Then, we discuss how to utilize the ergodicity of 
the Langevin process to show an almost Markovian behaviour of the spatial dynamics on long 
time scales. This eventually leads to a bound on the time scale on which the spatial dynamics 
is well approximated. 


4.1. Perturbation expansion 

The idea of perturbation expansion rests on the assumption that there exists a small problem 
parameter in which one can expand the objects of interest in a (formal) power series. As in 


Section 2.1 here this small parameter is the inverse of the damping in the Langevin dynamics, 


i.e. e := 7 1 where, for simplicity, we assume that the friction coefficient is scalar. For ease 
of presentation, we set the inverse temperature /3 = 1. 
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It turns out to be advantageous to work with the propagators (Koopman operators) instead 
of the transfer operators themselves. The difference is only of technical nature, since the prop¬ 
agators are the adjoints of the corresponding transfer operators. Denoting the propagators of 
the Langevin, Smoluchowski, and spatial dynamics by T^ an , Tg mol , and Tq, respectively, we 
have the explicit representations 


T Ln u(q,p) 

T Smo MQ) 

Tgw(q) 


^Lan Tjan\ 


E [u(q' t . -/p-f 


Lan _ 

% ~ <h 


E 


w(qf mol )\qr° l = Q 


Smol 


f E[u,(g t Lan ) 
Jv 


) I % Lan = Q, Po an = 


■n Lan 

P 0 


Lan 


= p\ 

p] fv{p)dp 


(41) 

(42) 

(43) 


where the expectation E[-] is taken with respect to the law of the stochastic forcing in the 
Langevin (for T/' an and Tq) and Smoluchowski (for Tg mol ) equations. The propagators T/' an 
and Tg mol are semigroups with generators 


Asmol = Ag-VgD-Vg 

A-Lan = P ‘ V q X7 qV • Vp + £ (Ap p • Vp) = AjJ am T £ Aou 

while Tq is not a semigroup, but ^Tq | = Ag mo i; in complete analogy with the theory 

presented above. 

To proceed, set A e := e~ 1 Ai, an . This scaling of ALan is called diffusive scaling and is due 
to the fact that the spatial dynamics gets slower and slower when friction is increased, and 
nontrivial dynamics only takes place on time scales of order e ^ 1 . The scaling of Ai^n by £ -1 


thus restores the relevant dynamics; see also (14). 


Now let (A e , u e ) be an eigenpair of A s , such that A e u £ = A £ u £ , and assume the existence of 
formal series expansions 


u £ = Uo + £U\ + £ 2 U 2 + • ■ • 

X e = Aq + £Ai + £^A 2 + ... 


It follows (see, e.g., [43] pp. 43], [39]) that uo(q,p) = uo(q), with Agmoi'«o = Ao^o, and u\(q,p) = 
p ■ X7gUo(q). This already gives a formal justification of the Smoluchowski dynamics as over¬ 
damped limit of the Langevin dynamics: on a time scale r = et (recall that Ai^n = eA e ) the 
position coordinate of the Langevin dynamics is governed by the Smoluchowski dynamics, up 
to fluctuations of order e. 


A closer look at the structure of the first terms in the eigenfunction expansion reveals even 
more. Metastability information is contained in eigenfuntions at nonzero eigenvalues, hence 
let Ao ¥= 0 7^ A e . Since Ag mo i is the formal adjoint of Lg mo i and AL an is the formal adjoint 
of -^Lan in C 2 , their eigenfunctions to different eigenvalues are orthogonal with respect to 
the corresponding scalar product. And, since the eigenfunctions of L^ an and Lsmol at the 
eigenvalue 0 are the canonical and Gibbs-Boltzmann densities / can and /q, respectively, we 
have that 



and 


fca.n(q,p)u 0 (q,p)dpdq = 

JqJv Jq 


fQ{q)uo{q)dq = o. 
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Being in the subspace orthogonal to / can , both functions decay exponentially under the action 
of the propagator. More precisely, let 


?y £ := max{ReA £ | 0 A A £ e a(A £ )} < 0 


be the real part of the nonzero eigenvalue of A e which is closest to zero; i.e. (£|r/ e |) _1 is 
the dominant time scale of the Langevin dynamics. Note that r] e = 0(1) as e —> 0, and 
lim sup e _, 0 Ve < 0- Now, both T^ &n u £ and Tp an 'Uo decay as ex.p(erj £ t) for t — > oo. We will 
utilize this with the perturbation expansion in the following computation. Its purpose is to 
estimate how far the Smoluchowski eigenfunction uq is from being an eigenfunction of the 
spatial propagator Tq. 


Tq u o(q) 


f ( T Lan«o) (q,p)fv(p)dp 
JV 

f { T Ln(.U e - (u £ -U 0 ))) (q,p)f V (p)dp 
JV 

e £Xet f u £ (q,p)f v (p)dp + 0 (e £1] A £ ) 

JV 

e £Xet f (u 0 (q) + eui(q,p) + 0(e 2 )) fv(p)dp + 0(e £Vst e) 
JV 

e £Xet uo(q) + 0(e eXet £ 2 ) + 0(e £rist E) as e —»■ 0 , 


where the third equality is obtained by utilizing the exponential decay of T Lni U e - U 0 ). The 
last equality follows from u\ and f-p being odd and even functions of p, respectively, hence 
the integral of their product vanishes. On the new, slower time scale t = et we obtain 

T £ ~ 1t u 0 = e XoT+ °^u 0 + e A0T+o(£) O(e 2 ) + e^ T 0(e). 

This means that uq is an approximate eigenfunction of the spatial propagator Tq 1t as long 
as e A ° r dominates the last two terms on the right hand side j£] It clearly dominates the second 
term (since we assume e to be small), hence we arrive at the desired condition by comparing 
it with the third0 

r < —— -r|l°g £ l or t <—— - r £ —1 1 log s\ (44) 

|Ao — 77e| |A 0 — Pel 

These estimates allow the following interpretation: 

1. While the standard result allows an approximation of the (position coordinate of the) 
Langevin dynamics by the Smoluchowski dynamics on a time scale et = t = 0(1) (as 
e —> 0), our estimate suggests that with respect to metastability analysis this time scale 
can be stretched by a factor | loge|. 


5 The spatial transfer operator is self-adjoint, hence normal. From the theory of pseudospectra for normal 
operators [JH] we know that if Tu = A u + ev for some linear operator T, u,v of modulus one, and A e R, 
then T has an eigenvalue in the e-neighborhood of A. 

6 Note that for x, y > 0 one has e~ x « e~ v if e y ~ x « 1, already achieved if y < x, meaning “y smaller than x 
up to some additive constant”. 
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2 . The more dominant an eigenvalue, i.e. the smaller |Ao — rje |> the longer the time scale 
is on which the Smoluchowski eigenmode approximates the corresponding eigenmode of 
the spatial propagator well. For the first subdominant eigenmode, where A e = rj £ , and 
hence Aq — % = 0(e), the estimate reads as r < e _1 | loge|, or equivalently, t < e -2 | loge|. 


In order to validate the estimate (44) also numerically, we perform the following experiment: 
Consider the Langevin system induced by the one-dimensional periodic potential 


V(q) = 1 + 3cos(27 rq) + 3cos 2 (27rq) — cos 3 (27rq) 


with constant mass matrix M = 1 at temperature (5 = 1. 

■ 10" 2 




Figure 1: The two wells of the periodic double well potential indicate two metastable regions 
in configuration space. The sign structure of the dominant eigenfunctions of T Q 
reveals them. 


For varying e = 7 -1 , we computed the largest lag time t = t v (e) such that the eigenfunctions 
u £ = u\ at the subdominant eigenvalue A £ = A l of Tq and Tg mol differ by less than a given 
threshold u, i.e. we compute 


L,(e) := inf < t > 0 


\\ u I( T q) - u l( T Smol)\\cj Q > " 


Figure [2] shows e >—> t u (e) for v = 0.05, and for comparison, the graph of e >—> c\ log(e)e 2 + C2 
(where we obtained the constants c\ and C2 by a least squares fit on the given data). Clearly, 


on the chosen domain for t u , there is an excellent agreement with the estimate (44). 


Although these first estimates allow merely a slight quantitative extension of the time scale 
on which the Smoluchowski dynamics approximates the spatial component of the Langevin 
dynamics well, it suggests that the consideration of further structural information from the 
perturbation expansion may allow for an extension of approximation time scales beyond the 
current, or inspire corrections terms to do so. 


4.2. The Mori-Zwanzig formalism 

In the previous section, we have analyzed a possibility to extend the time scales on which 
metastability information gained from the Smoluchowski equation is a good approximation 
to that of the actual model of interest, the Langevin dynamics. The argument, however, 
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Figure 2: e-dependence of the maximal lag time. The blue graph shows the largest lag time 
such that jj ul(Tg) — ul(Tg mo] )|| £ 2 < 0.05. The black graph is cilog(e)e~ 2 + C 2 

with ci % -1.04-10' 4 , C 2 « 1.07 • 10 _1 (from least-squares fitting). The eigenfunc¬ 
tions were computed using a simple Ulam discretization [2] of Tq and Tg mol with 
resolution 256. 


required the smallness of the inverse damping coefficient e = 7 -1 . In this section we turn to 
the question what can be done without this assumption. 

The Mori-Zwanzig representation decomposes the differential equation governing the state 
variable of interest into terms according to their dependence on the same quantities of inter¬ 
est. Note that the formalism itself is quite general [521 IS E]; we will give an brief introduce 
that is tailored to our needs; cf. also the related paper [45). 

Let A be the infinitesimal generator of some propagator semigroup, which we will formally 
denote by (e At )t^o- This propagator acts on scalar functions / which are functions of the full 

state x. Let x = ( x,x ), where x is the state of interest (also called the resolved variables). 

Let a distribution [i be given over the state space, and define the projection operator II as 
expectation with respect to // conditional on x: 

riff \ nff^ w [ff W A Ui x ) d T{x) 

n f{x) = n f(x) -.= [ f{x) | xj = ^ 

We are only interested in the evolution of average quantities conditional on x, i.e. in Iie tA . 
Note that in the conformational analysis setting x = (q,p), x = q, n is the canonical measure 
with density / can , A = ^Lan from above, and thus He tA = Tq, the spatial propagator. 

Let n 1 = Id — II denote the projection orthogonal on the space of functions of x. A 
modified Mori-Zwanzig representation yields 

— Ue tA = HAHe tA + UAU 1 f e^U 1 Ae^~ s)UA ds + LL iU L e mA . (45) 

dt Jo 

It can be obtained by applying Dyson’s formula nu on e tA in the second term on the right 
hand side of the identity ^He tA = nAIIe^ + HAH L e tA , where we tacitly assume that the 
orthogonal dynamics in the space of the unresolved variables is well-posed; see |18| for details. 
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Equation (45) is hard to interpret in this form. Again, the structure of the problem at hand 


aids us: the assertions of the following Lemma can be checked by direct computation. From 
now on we work in the conformational analysis setting, i.e. x, fi. and A are given as above. 

Lemma 5. Let f be a function independent of the variable p, i.e. f(q,p) = f(q) Vq,p. Then 
the following holds: 

(a) n/ = / and II A f = 0; 

(b) UAf = 0, thus e mA f = / Vt > 0, and UAU = 0; 

(c) UAU L Af = A Smol f. 


Applying identity (45) to some / being a function only of the spatial variable q, we obtain 
with Lemma 03 

,i rt 

ry = UAU 1 e sA U L Af ds . (46) 

at Jo 

Approximating the integral as ^ h(s)ds = th( 0) + 0(t 2 ), we get 

jTy = tUAU^Af + 0(t 2 ) = tAsmoif + 0(t 2 ). 


Integrating over t yields 


T y = f + -A Smol f + o(t 3 ). 


(47) 


This result can be seen as an analogue to (31), with the advantage that it has been derived from 


an identity , (46), which now offers the possibility of a systematic exploitation of quadrature 


rules of increasing order to approximate the integral on its right hand side. To this end, note 
that the range of II 1 is orthogonal to / can . Thus, by the considerations in the previous section, 
e tA U L u = Tf m U A u decays exponentially as t —> oo. This means that for some unknown A < 0 

ct As, 


the integral in (46) has the form e~ Xs g(s)e As ds , where g(s)e~ As = 0(1) as s —> oo. Integrals 

with exponential weights can be approximated by the Gaufi-Laguerre quadrature rule. 

There are two issues which have to be addressed: the unknown A < 0 and that we do not 
have an explicit expression for the function g(s ) = TL^II^tt for s > 0. For the former, one 
could use the eigenvalues of the Smoluchowski propagator as an initial guess, and try to re¬ 
fine this estimate subsequently. For the latter, bootstrapping techniques (just as used for the 
derivation of Runge-Kutta methods in numerical integration) could help us to approximate 
g(s) from derivatives of g at s = 0. 


We shall summarize which lines of attack so far the Mori-Zwanzig formalism offers to 
extend the time scales of approximation in molecular conformation analysis. 


1. Based on (46) we are able to state a second-order accurate approximation of the spatial 


propagator by the Smoluchowski propagator. Using higher order quadrature rules (either 
based on Taylor expansions of the integrand or Gaufi-Laguerre quadrature), it should be 
investigated whether simple higher order approximations can be derived as well. 
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2. Repeating the Mori-Zwanzig procedure starting from the decomposition 


jfIie tA = Ue tA UA 


Re tA U L A 


one arrives at different representations than (45). This is the usual approach (cf. [5]), as 
it allows the interpretation of the arising terms as “first order optimal prediction”, “mem¬ 
ory”, and “noise”, however does not lead as simply to (47) as (45) does. Nevertheless, it 
should be considered parallel to (45), as it may reveal other important characteristics of 
the spatial dynamics. 


3. In [5], the short memory approximation ^ h(s)ds « th( 0) has been used, however not for 
the position-momentum decomposition that we consider here. In the same work, different 
projections II have been considered, e.g. finite rank projections to a set of basis functions. 
The results of Section [471] suggest, that a projection of the spatial dynamics to the space 
spanned by dominant eigenfunctions of the Smoluchowski propagator may give a good 
approximation for larger times as well. 

4.3. Almost Markovian behaviour: on bounding the approximation time scales 

For small t , the non-Markovianity of spatial dynamics is an important feature which charac¬ 
terizes the density transport and metastability. We have seen that a A e ^(S 1 ) satisfies A —> 1 
as t —> 0 with a rate of 0(exp (— K,t 2 /2)) (for some suitable k), in contrast to the rate for 
semigroups of operators, which is 0(exp (—nt)). 

However, for larger t, S t exhibits a more regular, almost Markovian behaviour 13S3- We 
give ideas as to how this could be exploited for efficient computation of the eigenvalues of 5* 
in this time region. 

Relaxation times for momenta distributions. Langevin dynamics, the underlying model of 
S t , is both Markovian and ergodic [33]. Due to ergodicity, we observe the convergence of any 
density to the canonical density / can . Moreover, for sufficiently large damping, the relaxation 
of the momentum coordinates is significantly faster than of the position coordinates, which 
can be seen by considering the associated Fokker-Planck equation: 

3 f 1 

= (^Ham + If Lou) f , with L 0 U3 = J^v9 + • (g M ~ l p) ■ 

Thus, higher friction 7 implies that the Ornstein-Uhlenbeck-part dominates the time evolu¬ 
tion. The solution of the Ornstein-Uhlenbeck Fokker-Planck equation is 

g(t,p)= f K(t,p,r)g(0,r)dr, (48) 

Jv 

with 

K(t,p,r ) = (det(27r/3 _1 C(t))) exp ^ (p — re _7< ) T C'(t) -1 (p — re -7 *) 
and the covariance matrix 

C(t) = M (id — e _2 7 M_1 A . 
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Observe that the time variable t appears in K always multiplied by 7. Thus, the larger the 
damping 7, the more rapidly g(t, •) tends towards the stationary solution: 


lim K(t,p,r ) = f v (p) . 

t—> 00 

This suggests that we can find an optimal lag time r, such that for all t 0 r and for all 
/ : X -► R 

P Lj((hP) » /*0?)/can(<?,p) 

for some /* : Q —* M. 

We use this to argue in favor of “almost-Markovianity” of S t : In the following let t ^ r. 
For u : Q —> M there is an u t : Q —> M such that 


■^"Lan (■ u{q)fca,n{q,p )) « U t (q)f can {q,p). 


Using this and the semi-group property of P{ an , we get 

j PlL{ u (q)fcan(q,p))dp 

* ^PLn{ ut (q)fcan(q,p))dp 

= sV(?) 

* ^(/^y j p P Lan('w(^)/can(g,p))dp) 

= (S*) 2 «(g). 


Inductively, it follows that S ni « (S'*)” - for t 0 r, so in this sense, 5* is almost a semigroup 
for big enough t. As the relaxation rate in (48) scales with I/7, we expect the optimal lag 
time to do the same. 


Extrapolating the restored operator. Now assume that r > 0 is small enough to allow R T 
(or E T ) to be a reasonable approximation to S T . Then 


S (nr) ^ ( S Tjn ^ ( R ry 


We validate this with a simple numerical example. Using the one-dimensional periodic 


double-well potential introduced in Section |4. 1| , we want to compute the second largest eigen¬ 
value A 1 (S'*) < A °(S t ) = 1, which provides insights into the stability of the two metastable 
sets, as of Theorem [2] Note that A 1 (R t ) does not provide a good approximation for t signifi¬ 
cantly larger than r, as the error assymptotics of Corollary [2] only hold for t —> 0. 

With damping 7 = 5, a choice of r = I/7 = 0.2 seems reasonable, as by visual inspection, 
A 1 (S t ) in this region begins to show exponential decay. Moreover, R T and E T still approximate 
S T well enough: 


|A i (5 T ) - A (_R T )| « 0.15 , |A \S T ) - A \E T )\ « 0.12 . 


Figure 4.3 compares A 1 (S' t ) with A 1 (S T ) ,l , A 1 (i? T )” and A 1 (£’ T )" for n = 1,..., 10. As an 
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Figure 3: Dominant eigenvalue of S t and its approximations 


error estimate for the eigenvalues, we get 

\X(S nr ) - X((R T ) n ) | < |A(S nr ) - X({S T ) n ) | + |A((S T ) n ) - A((iT) n ) | 
^ |A(5 nr ) - A(5 r ) n | + |A(S T ) - A (R T ) \ n , 


with using the binomial formula to obtain the second inequality. The first term on the right 
hand side depends on the relaxation of the underlying process after lagtime r, and (for fixed 
n) decreases with increasing r. The second term depends on the approximation error of R T 
on S T and increases with increasing r. A balance between these two error sources must thus 
be found. Typically, the optimal lag time lies in the approximation region of R* and E l only 
for high damping 7 . This may or may not correspond to the physical model at hand, and is a 
significant limitation of the eigenvalue extrapolation method. Alternative restored operators 
(as proposed in Section 4.2) may allow the application for smaller values of 7 . 


5. Discussion 

We have considered the dynamics of the position coordinate for a molecular dynamics system 
given by the Langevin process in thermal equilibrium. After deriving the high friction limit 
in generalized coordinates, and obtaining the associated Kramers-Smoluchowski dynamics, 
we have seen that the Smoluchowski equations show up in the evolution of the position 
coordinates also for any 7 , in the short time asymptotics after rescaling time according to 
1 1 —> f 2 /2. This can be extended from position coordinates to essential and reaction coordinates 
(shown for the scalar case here). Finally, we discussed possible approaches to stretch the 
approximation time scales of these pseudo generator methods. Here we argued that these 
time scales on which a good approximation has to be achieved, are actually finite due to 
the ergodicity of the Langevin process, and their upper bound decreases with increasing 
damping 7 . 

The numerical experiments in [2] suggest that our theoretical findings on the asymptotic 
approximation error can be extended to the dominant spectrum as well, hence that the 
approach is applicable for metastability analysis. In order to be applicable to bio-chemically 
relevant systems, two main points have to be addressed: 
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Extension of the approximation quality for larger time scales. 

An important aspect of the pseudo generator approach is that it gives a practical tool 
to systematically derive coarse-grained models of molecular motion by projecting the 
dynamics onto a subspace of essential coordinates. It is yet unclear, however, how the 
projection onto essential coordinates influences the approximation quality of the pro¬ 
jected pseudo generator G| ss on the original process. We expect the dominant eigen¬ 
functions of G| ss to be usable to reliably identify metastable sets, if the selected essential 
coordinates are slow-moving in comparison to the ” non-essential” coordinates. A pertur¬ 
bation expansion in the style of Section 4.1 might be be able to provide a rigorous error 
estimate, and identify the role of the non-essential coordinates in the approximation. 
Also, the involvement of higher order derivatives of S* ss in the approximation scheme 
seems promising (cf. [3]). Taking into account higher order terms in the derivation of 
the pseudo generators seems especially useful when accurate coarse-grained diffusion 
models in terms of few collective variables are sought in cases when no explicit small 
parameter is available and therefore traditional averaging or homogenization methods 
to eliminate unresolved degrees cannot be applied [29j ,31]. 


(b) Numerical discretization. 

We derived a differential operator expression for projected pseudo-generators in essen¬ 
tial coordinates (G| ss , cf. Lemma[4]) and saw that this operator has a simple, closed form 
that can again be interpreted as the generator of a diffusion process. Its discretiza¬ 
tion, especially for multidimensional reaction coordinates, can be conveniently done via 
spectral collocation using the Feynman-Kac representation of the underlying partial 
differential equation or the transfer operator, depending on which type of problem is 
considered; for details, see e.g. ECU EH HU- We should stress, however, that while for 
the unprojected operator, G 2 , the collocation matrix can be set up analytically (2j, for 
the projected one, G| ss , high-dimensional integrals over non-essential degrees of freedom 
are involved. Sampling-based quadrature seems to be the natural treatment here (see 

mmm)- 

Further, even if the reduction to a comparatively small number of reaction coordinates 
can be carried out, the discretization of the corresponding pseudo generators will become 
computationally challenging due to the curse of dimension if there are more than, say, 
six of these degrees of freedom. On the other hand, the macroscopic dynamics of the 
molecular system is taking place on an essentially one-dimensional skeleton: Apart 
from motion within metastable states (i.e. the conformations of the molecule), the vast 
majority of conformational transitions occurs along a few dominant, low dimensional 
transition pathways mmm- While metastable states correspond to densities which 
are almost fixed points under the action of some transfer operator, the transitions can be 
modelled as curves in the space of densities. This picture alludes to numerical techniques 
for computing low-dimensional (invariant) sets in systems with higher dimensional state 
spaces 00, using ansatz functions of higher smoothness in combination with a meshfree 
approach M- 
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A. Coordinate expressions for the Smoluchowski equation 


In order to compute the right hand side in (12) explicitly, we observe that, for functions 

= {G^V) ■ Vlp 

where Vip(u, t ) is understood as the derivative with respect to the spatial argument, here u. 
Upon noting that 

-4ou (7"M = ~G~ l v , 

with Aq u acting component-wise from the left, we find that 


Anj^HamV’ = -(7 M ' V0 


- 1 . 


for the action of ^4 qu on ^HamV’ e ranTlou- Therefore 


^fHam^oU^Haml/’ — ^7^ 


M 


_d_ 

duj 


V +-vG~ 1 v 
_d_ 

dui v ' '•? 

d 2 ip 


dtp 

dui 


i,J L 

E (g'H Bu . ui . 


du* 


where upper indices indicate inverse matrices, i.e., y *- 7 = (7 1 )ij. Using the identity 

f v ■ Bv q u (v) dv = /j-tr {GB ), BeM. dxd , 

JRrf P 


with q u as given by (13), we can easily evaluate the integral in (12), which yields 

d 2 ip 




hi 


7 




dV 1 


duj 2/3 


-tdGW dip 1 (dyi d 


+ — tr G~ -— -A + - 


du* 


dui [3 \ duj du 


7 


y. 


dujUj / j 


In the last equation we have used the shorthand 

Alp = (^Ham^QuAHamVO Qui/v) dv , 

JV 

Employing Jacobi’s formula (det G)' = det G tv(G^ 1 G') and the fact that det G = det M det(V$ T V$), 
the above expression for A can be recast as desired: 

A = /3 _1 A — VU • V , 

where 


V = 7 and A = 


1 


.V 


det 7 7 X V ) , 


V / det 7 

denote gradient and Laplace-Beltrami operator with respect to 7 . Note that A no longer 
depends on the constant mass matrix M. 
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